In this example we will examine multivariate forecasting models using mvgam, which fits dynamic GAMs (DGAMs) using MCMC sampling (note that either JAGS or Stan is required; installation links are found here and here). First a simulation experiment to determine whether mvgam‘s inclusion of complexity penalisation works by reducing the number of un-needed dynamic factors. In any factor model, choosing the appropriate number of factors M can be difficult. The approach used by mvgam when sampling with JAGS is to estimate a penalty for each factor that squeezes the factor’s variance toward zero, effectively forcing the factor to evolve as a flat white noise process. By allowing each factor’s penalty to be estimated in an exponentially increasing manner (following Welty, Leah J., et al. Bayesian distributed lag models: estimating effects of particulate matter air pollution on daily mortality Biometrics 65.1 (2009): 282-291), we hope that we can guard against specifying too large a M. Note that when sampling with Stan, we capitalise on the superior sampling and exploration of Hamiltonian Monte Carlo to choose the number of factors by placing independent normal priors on factor standard deviations. Begin by simulating 6 series that evolve with a shared seasonal pattern and that depend on 2 latent random walk factors. Each series is 100 time steps long, with a seasonal frequency of 12. We give the trend moderate importance by setting trend_rel = 0.6 and we allow each series’ observation process to be drawn from slightly different Poisson distributions

set.seed(1111)
library(mvgam)
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-33. For overview type 'help("mgcv-package")'.
## Loading required package: parallel
dat <- sim_mvgam(T = 100, n_series = 6, n_lv = 2, 
    family = "poisson", mu_obs = runif(8, 
        4, 6), trend_rel = 0.6, train_prop = 0.85)
## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'tibble'
## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'pillar'
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

Have a look at the series

par(mfrow = c(3, 2))
for (i in 1:6) {
    plot(dat$data_train$y[which(as.numeric(dat$data_train$series) == 
        i)], type = "l", ylab = paste("Series", 
        i), xlab = "Time", bty = "L")
    box(bty = "L", lwd = 2)
}

par(mfrow = c(1, 1))

Clearly there are some correlations in the trends for these series. But how does a dynamic factor process allow us to potentially capture these dependencies? The below example demonstrates how. Essentially, a dynamic factor is an unobserved (latent) random process that induces correlations between time series via a set of factor loadings (\(\beta\)) while exercising dimension reduction. The loadings represent constant associations between the observed time series and the dynamic factor, but each series can still deviate from the factor through its error process and its associations with other factors (if we estimate >1 latent factor in our model).

A challenge with any factor model is the need to determine the number of factors M. Setting M too small prevents temporal dependencies from being adequately modelled, leading to poor convergence and difficulty estimating smooth parameters. By contrast, setting M too large leads to unnecessary computation. mvgam approaches this problem by formulating a prior distribution that enforces exponentially increasing penalties on the factor variances to allow any un-needed factors to evolve as flat lines. Now let’s fit a well-specified model for our simulated series in which we estimate random intercepts, a shared seasonal cyclic smooth and 2 latent dynamic factors

mod1 <- mvgam(data_train = dat$data_train, 
    data_test = dat$data_test, formula = y ~ 
        s(series, bs = "re") + s(season, 
            bs = c("cc"), k = 8), knots = list(season = c(0.5, 
        12.5)), use_lv = TRUE, n_lv = 2, 
    family = "poisson", trend_model = "RW", 
    use_stan = TRUE, chains = 4, burnin = 500)
## Info: Found int division at 'string', line 22, column 31 to column 48:
##   n_lv * (n_lv - 1) / 2
## Values will be rounded towards zero. If rounding is not desired you can write
## the division as
##   n_lv * (n_lv - 1) / 2.0

Look at a few plots. The estimated smooth function

plot_mvgam_smooth(object = mod1, series = 1, 
    smooth = "season")

And the true seasonal function in the simulation

plot(dat$global_seasonality[1:12], type = "l", 
    bty = "L", ylab = "True function", xlab = "season")

Check whether each factor was retained using the plot_mvgam_factors function. Here, each factor is tested against a null hypothesis of white noise by calculating the sum of the factor’s 1st derivatives. A factor that has a larger contribution to the series’ latent trends will have a larger sum, both because that factor’s absolute magnitudes will be larger (due to the weaker penalty on the factor’s precision) and because the factor will move around more. By normalising these estimated first derivative sums, it should be apparent whether some factors have been dropped from the model. Here we see that each factor is contributing to the series’ latent trends, and the plots show that neither has been forced to evolve as white noise

plot_mvgam_factors(mod1)

Now we fit the same model but assume that we no nothing about how many factors to use, so we specify the maximum allowed (the total number of series; 6). Note that this model is computationally more expensive so it will take longer to fit

mod2 <- mvgam(data_train = dat$data_train, 
    data_test = dat$data_test, formula = y ~ 
        s(series, bs = "re") + s(season, 
            bs = c("cc"), k = 8), knots = list(season = c(0.5, 
        12.5)), use_lv = TRUE, n_lv = 6, 
    family = "poisson", use_stan = TRUE, 
    trend_model = "RW", chains = 4, burnin = 500)
## If rounding is intended please use the integer division operator %/%.
## Info: Found int division at 'string', line 22, column 31 to column 48:
##   n_lv * (n_lv - 1) / 2
## Values will be rounded towards zero. If rounding is not desired you can write
## the division as
##   n_lv * (n_lv - 1) / 2.0

Use the same plots as model 1 to see if this model has also fit the data well

plot_mvgam_smooth(object = mod2, series = 1, 
    smooth = "season")

plot(dat$global_seasonality[1:12], type = "l", 
    bty = "L", ylab = "True function", xlab = "season")

Examining the factor contributions gives us some insight into whether we set n_lv larger than we perhaps needed to (with some of the factors clearly evolving as unconstrained, zero-centred random walks). These contributions can be interpreted similarly to ordination axes when deciding how many latent variables to specify

plot_mvgam_factors(mod2)

The very weak contributions by some of the factors are a result of the penalisation, which will become more important as the dimensionality of the data grows. Now onto an empirical example. Here we will access monthly search volume data from Google Trends, focusing on relative importances of search terms related to tick paralysis in Queensland, Australia

library(tidyr)
if (!require(gtrendsR)) {
    install.packages("gtrendsR")
}

terms = c("tick bite", "tick paralysis", 
    "dog tick", "paralysis tick dog")
trends <- gtrendsR::gtrends(terms, geo = "AU-QLD", 
    time = "all", onlyInterest = T)

Google Trends modified their algorithm for extracting search volume data in 2012, so we filter the series to only include observations after this point in time

gtest <- trends$interest_over_time %>% tidyr::spread(keyword, 
    hits) %>% dplyr::select(-geo, -time, 
    -gprop, -category) %>% dplyr::mutate(date = lubridate::ymd(date)) %>% 
    dplyr::mutate(year = lubridate::year(date)) %>% 
    dplyr::filter(year > 2012) %>% dplyr::select(-year)

Convert to an xts object and then to the required mvgam format, holding out the final 10% of observations as the test data

series <- xts::xts(x = gtest[, -1], order.by = gtest$date)
trends_data <- series_to_mvgam(series, freq = 12, 
    train_prop = 0.9)

Plot the series to see how similar their seasonal shapes are over time

plot(series, legend.loc = "topleft")

Now we will fit an mvgam model with shared seasonality and random intercepts per series. Our first attempt will ignore any temporal component in the residuals so that we can identidy which GAM predictor combination gives us the best fit, prior to investigating how to deal with any remaining autocorrelation. We assume a Poisson observation model for the response. Also note that any smooths using the random effects basis (s(series, bs = "re") below) are automatically re-parameterised to use the non-centred parameterisation that is necessary to help avoid common posterior degeneracies in hierarchical models. This parameterisation tends to work better for most ecological problems where the data for each group / context are not highly informative, but it is still probably worth investigating whether a centred or even a mix of centred / non-centred will give better computational performance. We suppress the global intercept as it is not needed and will lead to identifiability issues when estimating the series-specific random intercepts

trends_mod1 <- mvgam(data_train = trends_data$data_train, 
    data_test = trends_data$data_test, formula = y ~ 
        s(series, bs = "re") + s(season, 
            k = 8, m = 2, bs = "cc") - 1, 
    knots = list(season = c(0.5, 12.5)), 
    trend_model = "None", family = "poisson", 
    use_stan = TRUE, chains = 4, burnin = 1000)

Given that these series could potentially be following a hierarchical seasonality, we will also trial a slghtly more complex model with an extra smoothing term per series that allows its seasonal curve to deviate from the global seasonal smooth. Ignore the warning about repeated smooths, as this is not an issue for estimation.

trends_mod2 <- mvgam(data_train = trends_data$data_train, 
    data_test = trends_data$data_test, formula = y ~ 
        s(season, k = 8, m = 2, bs = "cc") + 
            s(season, series, k = 5, bs = "fs", 
                m = 1), knots = list(season = c(0.5, 
        12.5)), trend_model = "None", family = "poisson", 
    use_stan = TRUE, chains = 4, burnin = 1000)

How can we compare these models to ensure we choose one that performs well and provides useful inferences? Beyond posterior retrodictive and predictive checks, we can take advantage of the fact that mvgam fits an mgcv model to provide all the necessary penalty matrices, as well as to identify good initial values for smoothing parameters. Because we did not modify this model by adding a trend component (the only modification is that we estimated series-specific overdispersion parameters), we can still employ the usual mgcv model comparison routines

anova(trends_mod1$mgcv_model, trends_mod2$mgcv_model, 
    test = "LRT")
AIC(trends_mod1$mgcv_model, trends_mod2$mgcv_model)
summary(trends_mod1$mgcv_model)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## y ~ s(series, bs = "re") + s(season, k = 8, m = 2, bs = "cc") - 
##     1
## 
## Approximate significance of smooth terms:
##             edf Ref.df  Chi.sq p-value    
## s(series) 4.000      4   67947  <2e-16 ***
## s(season) 5.646      6 4464129  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.795   Deviance explained = 81.4%
## -REML = 1392.2  Scale est. = 1         n = 412
summary(trends_mod2$mgcv_model)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## y ~ s(season, k = 8, m = 2, bs = "cc") + s(season, series, k = 5, 
##     bs = "fs", m = 1)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.6755     0.4403   6.077 1.23e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                     edf Ref.df Chi.sq p-value    
## s(season)         5.507      6  177.7  <2e-16 ***
## s(season,series) 14.165     19 1851.9  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.831   Deviance explained =   85%
## -REML =   1322  Scale est. = 1         n = 412

Model 2 seems to fit better so far, suggesting that hierarchical seasonality gives better performance for these series. But a problem with both of the above models is that their forecast uncertainties will not increase into the future, which is not how time series forecasts should behave. Here we fit Model 2 again but specifying a time series model for the latent trends. We assume the dynamic trends can be represented using latent factors that each follow a RW process, and we will rely on the exponential penalties to help regularise any un-needed factors by setting n_lv = 4

trends_mod3 <- mvgam(data_train = trends_data$data_train, 
    data_test = trends_data$data_test, formula = y ~ 
        s(season, k = 8, m = 2, bs = "cc") + 
            s(season, series, k = 5, bs = "fs", 
                m = 1), knots = list(season = c(0.5, 
        12.5)), trend_model = "RW", use_lv = TRUE, 
    n_lv = 4, family = "poisson", use_stan = TRUE, 
    chains = 4, burnin = 1000)
## If rounding is intended please use the integer division operator %/%.
## Info: Found int division at 'string', line 24, column 31 to column 48:
##   n_lv * (n_lv - 1) / 2
## Values will be rounded towards zero. If rounding is not desired you can write
## the division as
##   n_lv * (n_lv - 1) / 2.0

Have a look at the returned Stan model file to see how the dynamic factors are incorporated

trends_mod3$model_file
##   [1] ""                                                                                                                     
##   [2] "//Stan model code generated by package mvgam"                                                                         
##   [3] "data {"                                                                                                               
##   [4] "int idx1[16]; // discontiguous index values"                                                                          
##   [5] "int idx2[4]; // discontiguous index values"                                                                           
##   [6] "int<lower=0> total_obs; // total number of observations"                                                              
##   [7] "int<lower=0> n; // number of timepoints per series"                                                                   
##   [8] "int<lower=0> n_lv; // number of dynamic factors"                                                                      
##   [9] "int<lower=0> n_sp; // number of smoothing parameters"                                                                 
##  [10] "int<lower=0> n_series; // number of series"                                                                           
##  [11] "int<lower=0> num_basis; // total number of basis coefficients"                                                        
##  [12] "real p_taus[1]; // prior precisions for parametric coefficients"                                                      
##  [13] "real p_coefs[1]; // prior locations for parametric coefficients"                                                      
##  [14] "vector[num_basis] zero; // prior locations for basis coefficients"                                                    
##  [15] "matrix[num_basis, total_obs] X; // transposed mgcv GAM design matrix"                                                 
##  [16] "int<lower=0> ytimes[n, n_series]; // time-ordered matrix (which col in X belongs to each [time, series] observation?)"
##  [17] "matrix[6,6] S1; // mgcv smooth penalty matrix S1"                                                                     
##  [18] "int<lower=0, upper=1> y_observed[n, n_series]; // indices of missing vs observed"                                     
##  [19] "int<lower=-1> y[n, n_series]; // time-ordered observations, with -1 indicating missing"                               
##  [20] "}"                                                                                                                    
##  [21] ""                                                                                                                     
##  [22] "transformed data{"                                                                                                    
##  [23] "// Number of non-zero lower triangular factor loadings"                                                               
##  [24] "// Ensures identifiability of the model - no rotation of factors"                                                     
##  [25] "int<lower=1> M;"                                                                                                      
##  [26] "M = n_lv * (n_series - n_lv) + n_lv * (n_lv - 1) / 2 + n_lv;"                                                         
##  [27] "}"                                                                                                                    
##  [28] ""                                                                                                                     
##  [29] "parameters {"                                                                                                         
##  [30] "// raw basis coefficients"                                                                                            
##  [31] "row_vector<lower=-30,upper=30>[num_basis] b_raw;"                                                                     
##  [32] ""                                                                                                                     
##  [33] "// dynamic factors"                                                                                                   
##  [34] "matrix[n, n_lv] LV;"                                                                                                  
##  [35] ""                                                                                                                     
##  [36] "// dynamic factor lower triangle loading coefficients"                                                                
##  [37] "vector[M] L;"                                                                                                         
##  [38] ""                                                                                                                     
##  [39] "// smoothing parameters"                                                                                              
##  [40] "vector<lower=0.0005>[n_sp] lambda;"                                                                                   
##  [41] "}"                                                                                                                    
##  [42] ""                                                                                                                     
##  [43] "transformed parameters {"                                                                                             
##  [44] "// GAM contribution to expectations (log scale)"                                                                      
##  [45] "vector[total_obs] eta;"                                                                                               
##  [46] ""                                                                                                                     
##  [47] "// trends and dynamic factor loading matrix"                                                                          
##  [48] "matrix[n, n_series] trend;"                                                                                           
##  [49] "matrix[n_series, n_lv] lv_coefs;"                                                                                     
##  [50] ""                                                                                                                     
##  [51] "// basis coefficients"                                                                                                
##  [52] "row_vector[num_basis] b;"                                                                                             
##  [53] ""                                                                                                                     
##  [54] "for (i in 1:num_basis) {"                                                                                             
##  [55] "b[i] = b_raw[i];"                                                                                                     
##  [56] "}"                                                                                                                    
##  [57] "// constraints allow identifiability of loadings"                                                                     
##  [58] "for (i in 1:(n_lv - 1)) {"                                                                                            
##  [59] "for (j in (i + 1):(n_lv)){"                                                                                           
##  [60] "lv_coefs[i, j] = 0;"                                                                                                  
##  [61] "}"                                                                                                                    
##  [62] "}"                                                                                                                    
##  [63] "{"                                                                                                                    
##  [64] "int index;"                                                                                                           
##  [65] "index = 0;"                                                                                                           
##  [66] "for (j in 1:n_lv) {"                                                                                                  
##  [67] "for (i in j:n_series) {"                                                                                              
##  [68] "index = index + 1;"                                                                                                   
##  [69] "lv_coefs[i, j] = L[index];"                                                                                           
##  [70] "}"                                                                                                                    
##  [71] "}"                                                                                                                    
##  [72] "}"                                                                                                                    
##  [73] ""                                                                                                                     
##  [74] "// derived latent trends"                                                                                             
##  [75] "for (i in 1:n){;"                                                                                                     
##  [76] "for (s in 1:n_series){"                                                                                               
##  [77] "trend[i, s] = dot_product(lv_coefs[s,], LV[i,]);"                                                                     
##  [78] "}"                                                                                                                    
##  [79] "}"                                                                                                                    
##  [80] ""                                                                                                                     
##  [81] "eta = to_vector(b * X);"                                                                                              
##  [82] "}"                                                                                                                    
##  [83] ""                                                                                                                     
##  [84] "model {"                                                                                                              
##  [85] "// parametric effect priors (regularised for identifiability)"                                                        
##  [86] "for (i in 1:1) {"                                                                                                     
##  [87] "b_raw[i] ~ normal(p_coefs[i], 1 / p_taus[i]);"                                                                        
##  [88] "}"                                                                                                                    
##  [89] ""                                                                                                                     
##  [90] "// prior for s(season)..."                                                                                            
##  [91] "b_raw[2:7] ~ multi_normal_prec(zero[2:7],S1[1:6,1:6] * lambda[1]);"                                                   
##  [92] ""                                                                                                                     
##  [93] "// prior for s(season,series)..."                                                                                     
##  [94] "for (i in idx1) { b[i] ~ normal(0, lambda[2]); }"                                                                     
##  [95] ""                                                                                                                     
##  [96] "for (i in idx2) { b[i] ~ normal(0, lambda[3]); }"                                                                     
##  [97] ""                                                                                                                     
##  [98] "// priors for smoothing parameters"                                                                                   
##  [99] "lambda ~ exponential(0.05);"                                                                                          
## [100] ""                                                                                                                     
## [101] "// priors for dynamic factor loading coefficients"                                                                    
## [102] "L ~ double_exponential(0, 1);"                                                                                        
## [103] ""                                                                                                                     
## [104] "// dynamic factor estimates"                                                                                          
## [105] "for (j in 1:n_lv) {"                                                                                                  
## [106] "LV[1, j] ~ normal(0, 1);"                                                                                             
## [107] "}"                                                                                                                    
## [108] ""                                                                                                                     
## [109] "for (j in 1:n_lv) {"                                                                                                  
## [110] "LV[2:n, j] ~ normal(LV[1:(n - 1), j], 1);"                                                                            
## [111] "}"                                                                                                                    
## [112] ""                                                                                                                     
## [113] "// likelihood functions"                                                                                              
## [114] "for (i in 1:n) {"                                                                                                     
## [115] "for (s in 1:n_series) {"                                                                                              
## [116] "if (y_observed[i, s])"                                                                                                
## [117] "y[i, s] ~ poisson_log(eta[ytimes[i, s]] + trend[i, s]);"                                                              
## [118] "}"                                                                                                                    
## [119] "}"                                                                                                                    
## [120] "}"                                                                                                                    
## [121] ""                                                                                                                     
## [122] "generated quantities {"                                                                                               
## [123] "vector[n_sp] rho;"                                                                                                    
## [124] "vector[n_lv] penalty;"                                                                                                
## [125] "matrix[n, n_series] ypred;"                                                                                           
## [126] "rho = log(lambda);"                                                                                                   
## [127] "penalty = rep_vector(1.0, n_lv);"                                                                                     
## [128] ""                                                                                                                     
## [129] "// posterior predictions"                                                                                             
## [130] "for(i in 1:n){"                                                                                                       
## [131] "for(s in 1:n_series){"                                                                                                
## [132] "ypred[i, s] = poisson_log_rng(eta[ytimes[i, s]] + trend[i, s]);"                                                      
## [133] "}"                                                                                                                    
## [134] "}"                                                                                                                    
## [135] "}"                                                                                                                    
## [136] ""

Inspection of the dynamic factors and their relative contributions indicates that the first factor is by far the most important

plot_mvgam_factors(trends_mod3)

Model 3 (with the dynamic trend) should provide far superior forecasts than relying only on the estimated smooths. Inspect the model summary

summary(trends_mod3)
## GAM formula:
## y ~ s(season, k = 8, m = 2, bs = "cc") + s(season, series, k = 5, 
##     bs = "fs", m = 1)
## 
## Family:
## Poisson
## 
## Link function:
## log
## 
## Trend model:
## RW
## 
## N latent factors:
## 4
## 
## N series:
## 4
## 
## N observations:
## 412
## 
## Status:
## Fitted using rstan::stan()
## 
## GAM smooth term estimated degrees of freedom:
##                     edf df
## s(season)         15.34  6
## s(season,series) 241.31  5
## 
## GAM coefficient (beta) estimates:
##                             2.5%          50%        97.5% Rhat n.eff
## (Intercept)          2.015451599  2.704062661  3.357046358    1  3827
## s(season).1         -0.725954600 -0.558499214 -0.391935085    1  4428
## s(season).2         -0.745968815 -0.530358930 -0.331730807    1  4622
## s(season).3         -0.454880908 -0.243435772 -0.034444696    1  4279
## s(season).4         -0.126528126  0.096704318  0.298432717    1  3892
## s(season).5          0.594950983  0.795557653  0.996021247    1  4303
## s(season).6          0.415311823  0.570308255  0.727358055    1  4434
## s(season,series).1  -0.025902484  0.112598184  0.270454054    1  4321
## s(season,series).2  -0.037057937  0.116948837  0.282879379    1  4343
## s(season,series).3  -0.235609594 -0.088417101  0.060802425    1  4144
## s(season,series).4   0.075293472  0.146413231  0.226387481    1  4321
## s(season,series).5  -0.246916278  0.481381022  1.224551401    1  2620
## s(season,series).6  -0.349067798 -0.141097228  0.039786395    1  4101
## s(season,series).7  -0.146865922  0.049291945  0.253029573    1  4909
## s(season,series).8  -0.197040972 -0.032679354  0.134856058    1  4284
## s(season,series).9  -0.090271168  0.001148285  0.094937042    1  4422
## s(season,series).10 -1.770867701 -0.993878954 -0.223043479    1  2947
## s(season,series).11 -0.086619972  0.065147423  0.223972260    1  4423
## s(season,series).12 -0.273598256 -0.104589717  0.054089139    1  4330
## s(season,series).13 -0.004182408  0.135707321  0.295333196    1  3979
## s(season,series).14  0.034888627  0.109534198  0.191690222    1  4289
## s(season,series).15 -0.929026883 -0.198378160  0.527048316    1  3431
## s(season,series).16 -0.308198606 -0.147532763 -0.002552626    1  4508
## s(season,series).17 -0.170926173 -0.010173041  0.139672936    1  4541
## s(season,series).18 -0.153108141 -0.006980581  0.139024325    1  4691
## s(season,series).19 -0.063331002  0.006514286  0.080009901    1  4444
## s(season,series).20  0.122248054  0.840516926  1.539518362    1  3712
## 
## GAM smoothing parameter (rho) estimates:
##                         2.5%          50%     97.5% Rhat n.eff
## s(season)          2.3903058  3.610948727  4.490441    1  4569
## s(season,series)  -2.5442110 -2.075032898 -1.560145    1  3100
## s(season,series)2 -0.7368928  0.006311793  1.167338    1  3516
## 
## [1] "n_eff / iter looks reasonable for all parameters"
## [1] "Rhat looks reasonable for all parameters"
## [1] "0 of 4000 iterations ended with a divergence (0%)"
## [1] "0 of 4000 iterations saturated the maximum tree depth of 12 (0%)"
## [1] "E-FMI indicated no pathological behavior"

Look at Dunn-Smyth residuals for some series from this preferred model to ensure that our dynamic factor process has captured most of the temporal dependencies in the observations

plot_mvgam_resids(trends_mod3, series = 1)

plot_mvgam_resids(trends_mod3, series = 2)

plot_mvgam_resids(trends_mod3, series = 3)

plot_mvgam_resids(trends_mod3, series = 4)

Perform posterior predictive checks to see if the model is able to simulate data that looks realistic and unbiased by examining simulated kernel densities for posterior predictions (yhat) compared to the density of the observations (y). This will be particularly useful for examining whether the Negative Binomial observation model is able to produce realistic looking simulations for each individual series.

ppc(trends_mod3, series = 1, type = "hist")

ppc(trends_mod3, series = 2, type = "hist")

ppc(trends_mod3, series = 3, type = "hist")

ppc(trends_mod3, series = 4, type = "hist")

Look at traceplots for the smoothing parameters (rho)

plot_mvgam_trace(object = trends_mod3, param = "rho")

Plot posterior predictive distributions for the training and testing periods for each series

plot_mvgam_fc(object = trends_mod3, series = 1, 
    data_test = trends_data$data_test)
## Out of sample DRPS:
## [1] 51.25285
## 

plot_mvgam_fc(object = trends_mod3, series = 2, 
    data_test = trends_data$data_test)
## Out of sample DRPS:
## [1] 13.42173
## 

plot_mvgam_fc(object = trends_mod3, series = 3, 
    data_test = trends_data$data_test)
## Out of sample DRPS:
## [1] 43.18961
## 

plot_mvgam_fc(object = trends_mod3, series = 4, 
    data_test = trends_data$data_test)
## Out of sample DRPS:
## [1] 33.33799
## 

Plot posterior distributions for the latent trend estimates, again for the training and testing periods

plot_mvgam_trend(object = trends_mod3, series = 1, 
    data_test = trends_data$data_test)

plot_mvgam_trend(object = trends_mod3, series = 2, 
    data_test = trends_data$data_test)

plot_mvgam_trend(object = trends_mod3, series = 3, 
    data_test = trends_data$data_test)

plot_mvgam_trend(object = trends_mod3, series = 4, 
    data_test = trends_data$data_test)

Given that we fit a model with hierarchical seasonality, the seasonal smooths are able to deviate from one another (though they share the same wiggliness and all deviate from a common ‘global’ seasonal function). Here we use the newdata argument to generate predictions for each of the hierarchical smooth functions (note that the intercept is still included in these plots so they do not center on zero)

newdat <- data.frame(season = seq(1, 12, 
    length.out = 100), series = levels(trends_data$data_train$series)[1])

plot_mvgam_smooth(object = trends_mod3, series = 1, 
    smooth = "season", newdata = newdat)

newdat <- data.frame(season = seq(1, 12, 
    length.out = 100), series = levels(trends_data$data_train$series)[2])

plot_mvgam_smooth(object = trends_mod3, series = 2, 
    smooth = "season", newdata = newdat)

newdat <- data.frame(season = seq(1, 12, 
    length.out = 100), series = levels(trends_data$data_train$series)[3])

plot_mvgam_smooth(object = trends_mod3, series = 3, 
    smooth = "season", newdata = newdat)

newdat <- data.frame(season = seq(1, 12, 
    length.out = 100), series = levels(trends_data$data_train$series)[4])

plot_mvgam_smooth(object = trends_mod3, series = 4, 
    smooth = "season", newdata = newdat)

Plot posterior mean estimates of latent trend correlations. These correlations are more useful than looking at latent factor loadings, for example to inspect ordinations. This is because the orders of the loadings (although constrained for identifiability purposes) can vary from chain to chain

correlations <- lv_correlations(object = trends_mod3)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
mean_correlations <- correlations$mean_correlations
mean_correlations[upper.tri(mean_correlations)] <- NA
mean_correlations <- data.frame(mean_correlations)
ggplot(mean_correlations %>% tibble::rownames_to_column("series1") %>% 
    tidyr::pivot_longer(-c(series1), names_to = "series2", 
        values_to = "Correlation"), aes(x = series1, 
    y = series2)) + geom_tile(aes(fill = Correlation)) + 
    scale_fill_gradient2(low = "darkred", 
        mid = "white", high = "darkblue", 
        midpoint = 0, breaks = seq(-1, 1, 
            length.out = 5), limits = c(-1, 
            1), name = "Trend\ncorrelation") + 
    labs(x = "", y = "") + theme_dark() + 
    theme(axis.text.x = element_text(angle = 45, 
        hjust = 1))

There is certainly some evidence of positive trend correlations for a few of these search terms, which is not surprising given how similar some of them are and how closely linked they should be to interest about tick paralysis in Queensland. Plot some STL decompositions of these series to see if these trends are noticeable in the data

plot(stl(ts(as.vector(series$`tick paralysis`), 
    frequency = 12), "periodic"))

plot(stl(ts(as.vector(series$`paralysis tick dog`), 
    frequency = 12), "periodic"))

plot(stl(ts(as.vector(series$`dog tick`), 
    frequency = 12), "periodic"))

plot(stl(ts(as.vector(series$`tick bite`), 
    frequency = 12), "periodic"))

Forecast period posterior predictive checks suggest that the model still has room for improvement:

ppc(trends_mod3, series = 1, type = "hist", 
    data_test = trends_data$data_test)

ppc(trends_mod3, series = 1, type = "mean", 
    data_test = trends_data$data_test)

ppc(trends_mod3, series = 2, type = "hist", 
    data_test = trends_data$data_test)

ppc(trends_mod3, series = 2, type = "mean", 
    data_test = trends_data$data_test)

ppc(trends_mod3, series = 3, type = "hist", 
    data_test = trends_data$data_test)

ppc(trends_mod3, series = 3, type = "mean", 
    data_test = trends_data$data_test)

ppc(trends_mod3, series = 4, type = "hist", 
    data_test = trends_data$data_test)

ppc(trends_mod3, series = 4, type = "mean", 
    data_test = trends_data$data_test)

Other next steps could involve devising a more goal-specific set of posterior predictive checks (see this paper by Gelman et al and relevant works by Betancourt for examples) and compare out of sample Discrete Rank Probability Scores for this model and other versions for the latent trends (i.e. AR2, AR3, Random Walk)